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Abstract — It  is  extremely  difficult  to  locate  the  con¬ 
tour  of  the  prostate  in  B-scan  ultrasound  images  auto¬ 
matically  by  computer  because  of  low  resolution  and  high 
noise  levels.  In  this  paper  we  present  a  semi-automatic 
prostate  contour  extration  scheme,  which  is  based  on  the 
wavelet  transform  and  active  contour  models,  or  snakes. 
The  ultrasound  image  is  first  decomposed  into  edge  maps 
at  different  resolutions  via  the  wavelet  transform.  Seed 
points  are  found  in  the  coarsest  edge  map  by  examining 
the  maxima  along  the  radial  profiles  which  emanate  from 
an  anchor  point  selected  manually.  These  seed  points  are 
used  to  initialize  a  snake,  which  will  evolve  across  the  edge 
maps  at  different  resolutions  and  eventually  converge  to 
the  contour  of  the  prostate. 
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I.  Introduction 

Prostate  cancer  has  become  a  focus  of  new  interest 
as  the  population  ages  and  the  death  rate  due  to  cardio¬ 
vascular  disease  decreases.  Advanced  prostate  cancer  is 
not  curable,  but  when  it  is  diagnosed  in  an  early  stage, 
it  is  curable  [1],  [2].  This  prompts  the  need  to  improve 
the  detection  rate  of  the  disease  in  early  diagnosis. 

Accurate  detection  of  prostate  contours  is  essential  to 
many  diagnostic  and  treatment  procedures.  For  example, 
prostate-specific  antigen  density  (PSAD)  is  very  useful  in 
differentiating  benign  from  malignant  prostate  disease. 
PSAD  involves  the  estimation  of  the  prostate  gland  vol¬ 
ume,  which  is  derived  from  the  boundaries.  Contour  in¬ 
formation  is  also  helpful  in  directing  biopsy  needles  to  the 
suspicious  site  and  monitoring  volume  change  of  hypo- 
echoic  lesions.  Currently  the  boundaries  of  the  prostate 
are  manually  outlined  by  a  urologist,  which  is  a  tedious 
and  challenging  task.  Moreover,  inter-observer  and  intra¬ 
observer  variability  is  inevitable.  So  automatic  contour 
extraction  by  computer  would  be  of  great  value. 

Unfortunately  it  is  extremely  difficult  to  determine 
the  prostate  boundaries  automatically  by  computer  due 
to  the  low  resolution  and  poor  contrast  of  B-scan  ultra¬ 
sound  images,  accompanied  by  the  high  level  of  speckle 
noise  due  to  scattering  and  other  complicated  interac¬ 
tions  between  ultrasonic  pulses  and  human  tissue.  More¬ 


over,  missing  boundary  segments  are  not  uncommon, 
which  are  attributed  to  acoustic  shadowing,  and/or  hy- 
perechoic  structures  within  or  around  the  prostate  gland. 

As  an  earlier  attempt,  a  contour  detection  scheme 
using  the  Laplacian  of  Gaussian  filter  was  applied  [3]. 
The  problem  with  this  scheme  was  that  it  is  difficult 
to  find  a  proper  threshold.  Too  many  or  too  few  edges 
may  be  found  depending  on  the  threshold  selected,  which 
makes  it  difficult  to  track  all  the  relevant  edges  that  cor¬ 
respond  to  the  object  and  form  a  closed  contour  from 
these  boundary  segments. 

Our  novel  contour  extraction  scheme  consists  of  the 
following  steps:  first  the  ultrasound  image  is  decomposed 
into  edge  maps  at  different  resolutions  using  the  wavelet 
transform.  Then  some  seed  points  are  found  in  the  coars¬ 
est  edge  map.  A  snake  is  initialised  using  these  seed 
points,  which  evolves  across  the  edge  maps  at  different 
scales  and  finally  converges  to  the  contour  of  the  prostate. 

II.  Contour  Extraction  Scheme 

A.  Multiscale  Edge  Detection 

Wavelets  are  families  of  scaled  and  shifted  versions 
of  a  mother  wavelet,  which  is  usually  a  waveform  of  lim¬ 
ited  duration  with  an  average  value  of  zero.  One  ma¬ 
jor  advantage  offered  by  wavelet  analysis  is  its  ability 
to  perform  local  analysis,  that  is,  the  wavelets’  localized 
support  means  that  some  local  features,  such  as  abrupt 
transitions  of  a  signal,  can  be  better  depicted  or  revealed 
by  wavelets  than  by  other  means. 

The  quadratic  spline  wavelet  and  its  scaling  func¬ 
tion  of  relevance  to  the  discrete  wavelet  transform  that 
follows  were  first  introduced  by  Mallat  and  Zhong  [4]. 
The  quadratic  spline  wavelet  is  the  first  order  deriva¬ 
tive  of  its  scaling  function,  which  closely  approximates 
a  Gaussian  function.  The  wavelet  decomposition  of  an 
image  I{x,y)  £  L2(R2)  at  scale  has  two  components, 
d^j{x:y)  and  d^(x,7/),  defined  by  the  convolution  with 
the  two  directional  wavelets  and  ^j(x,y)  re¬ 

spectively 

d%(x,y)  =  ip%(x,y)*I(x,y) 

d%j(x,y)  =  ip%(x,y)*I(x,y)  (1) 

The  two  directional  wavelets  are  given  by  the  partial 
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derivative  of  the  corresponding  spline  scaling  function 
02 as  follows 


d<hi  (x,  y) 

dx 

d<hi_  (x>  y ) 

dy 

By  substituting  (2)  into  (1)  we  get 

f  d%{x,y)  \  =  /  i>%{x,y)*I{x,y)  \ 

V  dXi  (x>  y)  )  V  ^ Xi  (x>  y)  * 1  y )  / 

*/)(*,  y)  \ 

*I)(x,y)  ) 

=  V(4)2j  *  I)(x,y) 


V’H  = 
= 


(2) 


(3) 


As  we  can  see  in  (3)  the  two  components  of  the  wavelet 
transform  of  an  image  I(x,y)  at  scale  are  equivalent  to 
the  gradient  of  the  image  smoothed  by  a  corresponding 
scaling  function  . 

The  modulus  of  the  wavelet  transform  is  given  by 


dZV (x,y)  =  \/\d%(x,  y)|2  +  \^ppy)\2  (4) 

and  the  angle  of  the  wavelet  transform  vector  is 

A2,  (*,»)=  (5) 

Shown  in  Fig.  1(a)  is  an  ultrasound  image  of  the  prostate, 
which  is  decomposed  up  to  scale  24  via  the  wavelet 
transform.  One  of  the  most  attractive  characteristics  of 
this  transform  is  that  the  notable  grey  scale  transitions 
within  the  original  image,  which  correspond  to  the  promi¬ 
nent  features  such  as  the  prostate  gland  boundaries,  are 
present  as  local  maxima  in  the  modulus  of  the  wavelet 
decompostition,  which  are  identified  as  bright  pixels  in 
Fig.l(b)-(c). 

The  transform  provides  us  with  the  edge  maps  of 
the  original  image  at  different  resolutions  from  the  coars¬ 
est  to  the  finest.  Note  the  contours  of  the  prostate  are 
prominent  in  the  coarser  (large  scale)  edge  maps.  In  con¬ 
trast,  the  noise  manifests  itself  in  a  random  pattern  and 
is  mostly  confined  to  the  edge  maps  of  small  scales.  This 
feature  implies  that  we  can  obtain  a  rough  boundary  of 
the  prostate  from  the  edge  map  at  the  coarsest  level  to 
avoid  the  interference  of  the  noise.  Then  we  can  refine 
the  boundary  by  incorporating  the  information  from  the 
edge  maps  of  lower  scales.  The  angle  of  the  two  compo¬ 
nents  of  the  wavelet  transform  also  contains  important 
information,  which  actaully  gives  us  the  orientation  of 
the  edge  vector  at  each  pixel.  The  angle  at  scale  24  is 
shown  in  Fig.  1(d).  The  information  of  both  the  modulus 
and  the  angle  will  be  employed  later  to  locate  the  seed 
points,  as  described  in  the  next  subsection. 


Fig.  1.  The  wavelet  decomposition  of  a  prostate  ultrasound  image, 
(a)  original  image,  (b)  modulus  (x,  y).  (c)  modulus  (x,  y). 
(d)  angle  A2 4  (x,  y) 

B.  Locating  the  Seed  Points 

The  next  step  is  to  locate  some  preliminary  boundary 
(seed)  points  in  the  coarsest  edge  map.  First,  an  anchor 
point  is  set  by  the  user,  who  also  indicates  the  rough 
position  of  the  first  seed  point.  A  number  of  radial  lines 
are  drawn  from  the  anchor  point  and  along  these  lines 
the  radial  profiles  are  extracted,  as  shown  in  Fig. 2 (a). 
Proceeding  radially  from  the  anchor  point,  a  seed  point 
is  always  characterised  by  a  local  maximum,  which  is 
larger  than  its  two  closest  neighbours  and  strictly  larger 


than  at  least  one  of  them. 

As  we  can  see  in  Fig. 2(b)  the  boundary  points  are 
present  as  the  local  maxima  along  the  radial  profiles  and 
the  loci  of  the  prostate  contour  are  vaguely  disernible. 
Most  of  the  desired  seed  points  are  of  large  amplitude. 
However,  some  of  them  are  very  weak  in  strength  re¬ 
sulting  in  broken  segments.  So  it  is  necessary  to  find  a 
strategy  which  somehow  favours  these  points  and  finally 
claims  them  as  seed  points  successfully.  As  we  can  see, 
every  two  adjacent  desired  seed  points  are  close  to  each 
other  because  of  the  fact  that  the  contour  is  relatively 
continuous.  Apparently  the  next  seed  point  should  be 
sought  within  a  small  neighbourhood  rather  than  “glob¬ 
ally”  along  the  radial  profile. 

Associated  with  each  edge  point  is  a  tangent  vector, 
which  is  perpendicular  to  the  orientation  of  the  edge  vec¬ 
tors,  that  is,  At2j{x^y)  =  A2j(x,y)  —  tt/2.  As  we  can  see 
in  Fig. 2(c),  the  tangent  vector  of  the  previous  seed  point 
is  roughly  pointing  to  the  position  of  next  desired  seed 
point.  We  also  take  the  perpendicular  line  of  the  next 
radial  line.  The  intersection  points  of  these  two  lines 
with  the  next  radial  line  defines  a  neighbourhood  within 
which  the  next  desired  seed  point  will  be  sought.  The 
points  outside  the  neighbourhood  are  adversely  weighted 
so  that  they  will  not  challenge  the  local  maximum  within 
the  neighbourhood,  which  is  then  asserted  as  the  next 
seed  point,  as  shown  in  Fig. 2(d). 

Fig. 3  shows  all  the  seed  points  successfully  identified 
using  this  strategy.  The  plus  sign  indicates  the  anchor 
point  and  the  cross  is  the  first  edge  point,  both  selected 
by  the  user.  These  seed  points  will  be  used  to  initialise 
the  snakes  which  are  expected  to  find  the  contour  even¬ 
tually. 

C.  Snakes 

A  snake  is  a  dynamic  curve, 

v(s)  =  (x(s),y(s)),  s  £  [0, 1] 

which  moves  within  an  image  trying  to  minimize  its  en¬ 
ergy,  which  is  defined  as  the  weighted  sum  of  an  internal 
energy  and  an  external  energy  [5],  [6]. 

E  snake  =  j  intip)  T  If-^ext  {v^)ds 

Jo 

=  jA(a\v'(8)\2+f}\v"{8)\*) 

+  7  Eext{y)ds  (6) 

The  internal  energy  comes  from  within  the  curve  itself 
and  is  defined  as  energy  of  continuity  and  energy  of  curva¬ 
ture.  The  external  energy  is  usually  set  as  the  gradient  of 
the  image,  which  will  draw  the  snake  towards  the  salient 
features  of  the  image,  such  as  edges  (high  amplitude  of 
gradient).  The  behaviour  of  the  snake  is  controlled  by 


the  regularization  parameters  a,/?, 7.  A  large  a  discour¬ 
ages  stretching  and  makes  the  snake  behave  like  an  elastic 
string,  which  can  shrink  or  contract.  A  snake  with  large 
/?,  which  discourages  bending  and  gives  a  dominant  effect 
of  smoothing,  behaves  more  like  a  rigid  rod.  If  we  set  f3 
to  zero,  that  is,  switch  off  the  influence  of  the  bending 
effect,  we  can  allow  a  sharp  corner  to  form.  The  energy 
minimizing  process  can  also  be  interpreted  from  a  point 
of  view  of  force  balance.  When  the  elastic  force  (cor¬ 
responding  to  the  energy  of  continuity),  bending  force 
(corresponding  to  the  energy  of  curvature),  and  poten¬ 
tial  force  (corresponding  to  the  external  energy)  reach  an 
equilibrium,  the  snake  is  expected  to  rest  on  the  objects 
of  interest  within  the  image. 

In  our  contour  extraction  scheme,  the  external  energy 
is  designated  as  the  modulus  of  the  wavelet  transform, 
that  is,  Eext(v )  =  d2jv(v).  The  seed  points  obtained 
from  the  previous  step  are  used  to  initialise  a  snake  which 
evolves  across  the  edge  maps  at  different  resolutions  and 
finally  stabilises  and  rests  on  the  contour  of  the  prostate, 
as  shown  in  the  next  section. 

III.  Results 

The  seed  points  identified  at  scale  24  (Fig. 3)  are  used 
to  intialise  a  snake,  which  starts  evolving  from  the  edge 
map  at  scale  23  (Fig.  1(b)).  The  control  parameters  are 
a  =  0.2,  f3  =  0.2.  7  is  set  to  1.1  to  emphasize  the  im¬ 
age  force  at  higher  scales  and  reduced  to  1.0  to  make  the 
snake  immune  to  the  interference  from  the  noise.  As  we 
can  see  in  Fig. 4,  the  snake  first  contracts  and  smooths  it¬ 
self  due  to  the  influence  of  its  internal  energy.  Gradually 
it  is  attracted  to  the  prominent  edges  within  the  image 
due  to  the  potential  force.  The  stablised  snake  at  scale  23 
is  then  used  to  intialise  a  snake  which  will  evolve  within 
the  edge  map  at  scale  2 2  and  directs  itself  towards  finer 
feature  within  the  edge  map.  The  final  contour  found  by 
the  snakes  is  shown  in  Fig. 5,  which  closely  matches  the 
edges  disernible  in  the  original  image. 
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Fig.  2.  Locating  the  seed  points  by  maxima  in  radial  profiles,  (a) 
emanating  radial  lines,  (b)  radial  profiles,  (c)  the  tangent  vectors 
and  the  perpendiculars,  (d)  locating  the  maximum. 


Fig.  4.  Snake  evolving  at  scale  23 


Fig.  5.  Contour  extracted 


